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Abstract.  In  this  note,  we  present  a  method  for  flattening  anatomi¬ 
cal  surfaces  such  as  branched  vessels  and  intestinal  tracts  in  an  area¬ 
preserving  way.  This  method  is  based  on  the  theory  of  optimal  mass 
transport  and  conformal  mapping  of  surfaces.  The  flattened  representa¬ 
tions  differ  minimally  from  conformality  in  a  certain  precise  sense.  Po¬ 
tential  applications  include  the  detection  and  visualization  of  pathologies 
such  as  stenoses  and  polyps. 


1  Introduction 

Surface  deformations,  in  particular  the  flattening  of  highly  undulated  and 
branched  surfaces,  is  an  active  area  of  research  in  the  field  of  medical  imag¬ 
ing  and  visualization,  For  example,  flattened  representations  of  brain  surface  are 
important  in  the  study  of  neural  activities  within  the  three  dimensional  folds  of 
brain  surface;  see  [HP]  and  the  references  therein.  Flattened  versions  of  the  colon 
surface  may  be  a  good  complement  for  CT  colonography  or  virtual  colonoscopy; 
see  mm  and  the  references  therein. 

There  have  been  a  number  of  approaches  to  surface  flattening.  Tn[1  2|  ,  a  vi¬ 
sualization  technique  is  proposed  which  uses  cylindrical  and  planar  map  pro¬ 
jections.  Methods  based  on  quasi-isometric  and  quasi-conformal  flattenings  of 
brain  surface  have  been  considered  in  m,  among  others.  Wang  m  presented 
a  method  for  unravelling  colon  surface  based  on  electrical  field.  Typically,  these 
methods  do  not  guarantee  bijectivity  of  the  mapping.  Algorithms  based  on  con¬ 
formal  mapping  have  also  been  applied  on  this  problem  mm-  These  methods 
preserve  angles,  and  in  this  sense  preserve  local  geometry.  They  also  guarantee 
bijectivity.  However,  these  conformal  flattening  methods  are  not  area-preserving, 
since  a  mapping  cannot  be  both  angle-preserving  and  area-preserving  unless  the 
original  surface  has  zero  Gaussian  curvature.  This  problem  becomes  particularly 
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pronounced  when  we  wish  to  construct  a  flattened  representation  for  a  multi- 
branched  surface. 

Here  we  take  another  approach  based  on  the  theory  of  optimal  mass  trans¬ 
port.  Our  method  makes  minimal  adjustments  to  an  initial  conformal  mapping 
to  find  an  area-preserving  mapping.  [Tj. 

2  Conformal  Flattening  of  a  Multi-branched  Vessel 

We  now  briefly  review  our  approach  for  constructing  a  conformal  flattening  of 
a  multi-branched  vessel.  The  algorithm  contains  two  steps:  dividing  the  whole 
vessel  into  several  Y-shaped  segments  using  the  harmonic  skeleton,  and  then 
flattening  each  Y-shaped  segment.  More  precisely,  our  method  is  based  on  the 
solution  of  a  harmonic  equation  on  the  surface,  from  which  a  “harmonic  skele¬ 
ton”  or  medial  axis  is  derived.  The  harmonic  skeleton  is  easy  to  calculate,  is 
guaranteed  to  be  smooth,  and  can  be  used  as  a  “center  line”  for  a  fly-through. 

Let  £  be  a  triangulated  surface,  without  self-intersections,  which  is  topolog¬ 
ically  a  tube  with  several  branches.  A  branched  blood  vessel  is  an  example  of 
such  a  surface.  We  assume  each  branch  terminates  in  a  boundary  loop  which  we 
call  Oi  (i  =  0 ...TV).  The  first  step  in  our  algorithm  is  to  solve  a  Dirichlet  prob¬ 
lem  Au  =  0  on  the  interior  of  S  with  proper  boundary  conditions.  One  of  the 
boundaries,  (Jq,  is  selected  to  be  the  root.  The  value  of  u  on  the  root  boundary  is 
assigned  to  be  0.  Given  this  0  contour,  there  exists  a  set  of  triangles  whose  base 
vertices  are  coincident  with  the  0  (starting)  contour  and  whose  apex  vertices 
describe  the  successive  contour.  This  description  can  be  repeated  for  n  pairs  of 
triangle  sets  until  the  successive  contour  reaches  a  boundary  other  than  root. 
The  value  of  u  on  this  boundary  is  then  assigned  to  be  n.  The  equation  Au  =  0 
is  solved  using  standard  finite  element  techniques  [7]. 

The  second  step  is  to  build  a  tree-like  structure  we  call  the  harmonic  skeleton. 
Using  the  function  u  solved  for  in  the  first  step,  we  find  level  sets  on  the  surface, 
i.e.  sets  {x\ u(x)  =  is}  for  values  of  is  ranging  from  0  to  the  maximum  of  u.  We 
partition  each  level  set  into  classes  according  to  connectivity  on  the  surface.  The 
centroid  of  each  class  of  points  corresponds  to  a  point  on  the  harmonic  skeleton. 
By  increasing  is  from  0  to  its  maximum,  we  can  build  a  structured  tree.  The 
locations  of  the  bifurcations  along  the  tree  are  obvious  since  it  is  easy  to  see 
where  one  centroid  splits  into  two. 

The  third  step  is  to  refine  the  harmonic  skeleton  by  replacing  the  boundary 
values  used  previously  in  step  1  with  new  values  obtained  as  the  distance  along 
the  skeleton  from  the  root  boundary  to  the  other  boundaries.  We  again  solve 
the  corresponding  Dirichlet  problem  and  extract  a  refined  harmonic  skeleton.  It 
is  easy  to  cut  a  branched  vessel  into  several  Y-shaped  segments  with  the  help  of 
this  construction. 

Details  for  the  conformal  flattening  of  a  Y-shaped  vessel  (i.e.,  with  only  a 
single  branch  point)  have  been  presented  in  |16|  to  which  we  refer  the  interested 
reader.  Essentially,  the  method  involves  solving  another  Dirichlet  problem  to 
obtain  a  harmonic  conjugate  for  u. 
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It  is  well  known  that  conformal  flattenings  are  “similarities  in  the  small”  in 
the  sense  that  over  small  regions,  the  flattened  version  differs  from  the  original 
by  a  scale  factor  alone.  As  a  result,  the  flattened  surface  retains  much  of  the 
visual  quality  of  the  original.  However,  conformal  mappings  do  not  in  general 
preserve  surface  area.  Indeed,  some  areas  of  the  surface  may  be  greatly  enlarged 
or  reduced.  The  problem  becomes  exacerbated  when  there  are  a  large  number  of 
branches.  To  address  this  problem,  we  consider  adjusting  the  conformal  mapping 
in  a  minimal  way  so  as  to  correct  the  distortion  of  area.  By  making  only  these 
minimal  adjustments,  we  hope  to  preserve  the  desirable  qualities  of  the  conformal 
mapping  as  much  as  possible. 

3  Introduction  Optimal  Mass  Transport 

3.1  Background 

Our  method  for  minimally  adjusting  a  conformal  mapping  is  based  on  the  theory 
of  optimal  mass  transport,  also  referred  to  as  the  Monge-Kantorovich  problem. 
Here  we  present  a  review  of  the  basic  theory.  Let  l?o  and  f2 1  be  two  domains 
in  R2,  having  smooth  boundaries.  On  each  of  these  domains  we  assume  we  are 
given  a  priori  a  mass  density  function  p q  and  Hi  respectively.  We  assume  that 
the  same  total  mass  is  associated  with  each  of  these  densities. 

We  will  be  considering  a  class  of  diffeomorphisms  u  from  f?0  to  fl\  which 
satisfy  the  equation 

Ho  =  \Du\pi  o  u.  (1) 

Here  Du  is  the  matrix  of  first  derivatives  of  u ,  and  o  represents  composition  of 
functions.  This  is  the  “Jacobian  equation”,  which  is  a  mass  preservation  (MP) 
constraint.  For  a  function  satisfying  (TJ)  we  will  write  u  £  MP.  The  special  case 
where  /n  =  1  is  of  particular  interest  for  the  current  application.  In  this  case 
Equation  CD  reduces  to  Ho  =  \Du\,  and  we  see  that  this  equation  is  simply  a 
hard  constraint  on  the  Jacobian  \Du\  of  u,  the  Jacobian  giving  the  amount  of 
scaling  of  area  that  occurs  locally  under  the  mapping  u. 

We  are  interested  in  finding  an  MP  mapping  u  which  differs  minimally  from 
the  identity.  To  this  end,  we  introduce  a  squared  L2  Kantorovich- Wasserstein 
penalty  functional  on  u  £  MP,  defined  as: 

M[u)  :=  [  ||u(a;)  —  x\\2  Ho(x)dx  (2) 

J  £2q 

This  functional  m  is  seen  to  place  a  penalty  on  the  distance  the  map  u  moves 
each  bit  of  material,  weighted  by  the  material’s  mass.  Hence,  the  Kantorovich- 
Wasserstein  metric  may  be  thought  of  as  the  cost  associated  with  transporting, 
via  u,  a  distribution  of  material  given  by  Ho-  The  resulting  distribution  of  ma¬ 
terial  is  constrained  to  be  the  given  density  fi\ ,  according  to  ([T]). 

Theoretical  results  EEI  show  that  there  is  a  unique  minimizer  u  £  MP  of 
M,  and  that  this  minimizer  is  characterized  as  being  the  gradient  of  a  convex 
function  w,  i.e.,  u  =  Vie. 
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3.2  Finding  the  Minimizer 

There  have  been  many  approaches  for  finding  the  solution  u  to  the  optimal  mass 
transport  problem.  In  |ll8|.  a  new  approach  for  finding  the  optimal  map  was 
proposed,  based  on  the  notion  of  polar  factorization  |6I3|.  Here  we  present  the 
main  idea  behind  this  work. 

Our  algorithm  for  finding  u  consists  of  two  steps.  The  first  step  is  to  find 
an  initial  mapping  u  from  to  1? i  satisfying  the  Jacobian  equation  0  .  The 
second  step  is  to  adjust  this  mapping  iteratively  using  gradient  descent  in  or¬ 
der  to  minimize  the  functional  M,  while  constraining  u  so  that  it  continues  to 
satisfy  QD-  a  method  for  finding  an  initial  mapping  is  described  below.  Here, 
we  derive  the  appropriate  gradient  descent  equation,  considering  a  more  general 
cost  functional  M  given  by 

M[u4]  =  f  d>{ut( x )  —  x)p,o(x)dx,  (3) 

J  J?0 

where  <P  :  R2  ->  R  is  a  positive  C 1  function.  The  L2  Kantorovich- Wasserstein 
functional  m  corresponds  to  the  case  <I>(x)  =  ||x||2. 

We  think  of  u  as  a  function  of  time  t,  and  introduce  the  notation  u4,  with  u° 
being  the  initial  mapping.  Each  such  u 4  can  be  written  as  the  composition 

w4  =  u°  o  (s4)-1,  or  (4) 

u°  =  utos\  (5) 

where  s4  is  a  mass  preserving  mapping  from  f?o  to  itself  satisfying  /io  =  |  Dst\  po° 
s4.  The  derivative  of  s4  with  respect  to  time  is  described  by  a  vector  field  v 4  on 
f?o,  i.e. 


o  s 


(6) 


where  the  vector  field  v 4  must  satisfy  div(/ro  u4)  =  0  in  order  to  preserve  the 
density  /Lto-  By  the  chain  rule  applied  to  Equation  (ED.  the  evolution  of  u4  must 
satisfy 

^L+iVV  =  0.  (7) 

By  a  change  of  variables  x  =  s4(y)  in  ([3D,  together  with  (0D,  we  find 


M(t)=  [  ${u°{y)  -  s^y))  p0(y)dy.  (8) 

Jn 

From  this  one  computes  that  the  cost  decreases  according  to 
rlM  t  f 

—  =  -  J  (V0(u°  -  s),  —)yody  =  -  J  (W(w4(a:)  -  x),ii0vt)  dx. 

Clearly,  were  it  not  for  the  constraint  div(p,o  v4)  =  0,  we  could  take  /io  c4  = 
\/d>(ut(x)—x)  to  decrease  M.  Instead,  we  take  yo  vt  to  be  the  divergence-free  part 
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of  V^(it*(x)  —  x ),  finding  this  part  through  the  use  of  Helmholtz  decomposition. 
This  step  involves  solving  Laplace’s  equation  with  boundary  conditions  set  to 
keep  the  flow  of  confined  in  f?0-  The  resulting  equation  for  updating  ul  has 
the  following  form: 

f)  o/t  1 

—  = - Du1  ■  (I  -  VZ\-1div)W(w*  -  id),  (9) 

at  po 

where  id  denotes  the  identity  map,  and  I  is  the  identity  matrix.  For  more  details 
on  this  technique,  please  refer  to  EHH3- 


3.3  Finding  an  Initial  MP  Mapping 

A  general  method  for  finding  an  initial  mapping  for  irregularly  shaped  domains 
can  be  found  in  the  work  of  Moser  m-  In  our  algorithm,  both  source  and  target 
domains  are  rectangles  of  the  same  size.  Accordingly,  we  can  use  the  following 
simple  method,  as  described  in  m- 

Assume  l?o  =  [0,  A]  x  [0,  B]  and  Q\  =  [0,  A]  x  [0,  B],  Our  initial  mapping  will 
be  of  the  form  u°(x,y)  =  (a(x),  b(x,  y)).  Since  both  po  and  Hi  are  positive,  we 
may  define  a(x )  implicitly  by  the  following  equation: 

pa(x)  pB  px  pB 

/  /  MiO?  ,y)dydr]=  /  /  p0(ri,y)dydr].  (10) 

J o  Jo  Jo  Jo 


Note  that  a(0)  =  0  and  a(A)  =  A  by  the  assumption  that  both  domains  contain 
same  amount  of  mass.  By  taking  the  derivative  respect  to  x,  we  find 

Mi  (a(x),y)dy=  /  /Jj0(x,y)dy.  (11) 

Jo 

We  may  now  therefore  define  b  =  b{x,y)  implicitly  by  the  equation 


rb(x,y)  py 

a'(x )  /  Hi(a(x),p)dp=  /  po(x,p)dp.  (12) 

Jo  Jo 

Taking  the  derivative  with  respect  to  y  shows  that  u°  satisfies  the  constraint  ©. 
This  construction  of  u°  can  be  interpreted  as  finding  the  optimal  mass  transport 
in  the  x  direction,  then  in  the  y  direction  for  each  x. 


4  Optimal  Transport  Applied  to  Conformal  Mapping 

Our  goal  is  to  apply  the  optimal  transport  theory  to  adjust  a  conformal  mapping 
minimally  in  order  to  obtain  an  area  preserving  mapping.  Let  us  now  assume  we 
have  constructed  a  conformal  mapping  /  of  a  triangulated  surface  to  a  region 
range (/)  in  the  plane. 

For  simplicity,  and  without  loss  of  generality,  we  assume  the  initial  surface 
and  rang e(/)  have  areas  equal  to  1.  We  can  always  scale  homothetically  to  make 
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this  so.  Define  /t o  on  range  (/)  to  be  the  area  of  a  triangle  on  the  original  surface 
divided  by  the  the  area  of  the  triangle  once  flattened.  Extend  the  function  /t q 
to  a  rectangular  region  170  surrounding  this  range  by  setting  /r0  to  1  outside  of 
range(/).  Figure [4] shows  such  a  function  /to,  the  dark  color  representing  enlarged 
areas  and  the  light  color  representing  reduced  areas.  Let  fl\  =  J70,  and  set  /ti  to 
be  uniformly  1  throughout  !?i. 

We  now  solve  the  optimal  transport  problem  as  described  in  the  previous 
section.  Since  /ti  =  1,  the  constraint  0  reduces  to  /To  =  \Du\,  and  so  by  the 
construction  of  /To  a  mapping  u  which  satisfies  this  constraint  compensates  ex¬ 
actly  for  the  distortion  in  area  that  occurred  during  flattening.  Further,  since 
we  minimize  the  functional  to  find  the  minimizing  u,  our  solution  differs 
minimally  from  the  identity. 

5  Implementation  and  Examples 

We  use  standard  techniques  to  solve  Equation  ©.  In  particular  we  have  em¬ 
ployed  an  upwinding  scheme  when  computing  Du*  7  and  the  FFT  when  inverting 
the  Laplacian  on  a  rectangular  grid.  Standard  centered  differences  were  used  for 
other  spatial  derivatives.  Once  we  numerically  solve  for  the  right  hand  side  of 
we  use  the  result  to  update  u* .  The  optimal  map  is  obtained  as  t  — >  oo.  In 
practice,  we  iterate  until  convergence  with  respect  to  a  specified  tolerance. 

Our  first  example  is  a  dataset  provided  by  the  Surgical  Planning  Lab  of 
Brigham  and  Women’s  Hospital,  Harvard  Medical  School.  The  dataset  is  a  vessel 
surface  extracted  from  a  256  x  256  x  47  MRA  brain  image.  By  using  a  segmenta¬ 
tion  method  based  on  active  contour  model,  we  found  the  surface  of  a  section  of 
vessel.  We  then  generated  a  triangulation  of  the  surface  using  the  Visualization 
Toolkit  as  shown  in  Figure  [1]  and  triangles  at  the  ends  were  removed  to  create  a 
tubular  surface.  Next,  we  used  the  technique  described  in  Section^  we  produced 
a  conformal  flattening  of  this  surface  (Figure  [2]),  together  with  a  harmonic  skele¬ 
ton  (Figure  El).  Then,  as  described  in  Section  [I]  density  maps  were  generated 
(Figure  EJ  based  on  the  distortion  of  area  due  to  flattening.  Figure  E]  presents  a 
histogram  of  this  distortion,  from  which  it  can  be  seen  that  many  triangles  were 
enlarged  or  reduced.  Finally,  an  area-preserving  mapping  (Figure  ED  was  con¬ 
structed  based  on  the  algorithm  described  in  Section  EOl  The  computation  took 
about  9  minutes  for  the  conformal  flattening  and  about  20  minutes  for  the  area 
preserve  flattening,  both  on  a  1  GHz  Linux  system.  All  surfaces  in  the  figures 
are  shaded  by  using  outward  normal  vectors  from  the  original  surface. 

The  second  example  is  a  dataset  provided  by  Emory  Hospital.  The  dataset 
is  a  512  x  512  x  91  CT  image  for  carotid  vessel.  The  flattening  approach  for  this 
dataset  was  similar  to  that  of  the  first  one,  except  that  there  was  only  one  branch 
point  and  hence  it  was  not  necessary  to  generate  a  harmonic  skeleton.  Figure  |7] 
is  the  segmented  surface  and  Figure  [H] is  the  area-preserving  presentation. 

Color  images  of  these  examples  can  be  found  at: 
http://www.prism.gatech.edu/~gte538w/miccai2003/miccai2003.htm 
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Fig.  1.  A  segmented  brain  vessel 


Fig.  5.  The  statistics  for  Fig  [31 


Fig.  6.  The  area-preserving  mapping  of 
Fig.  2.  The  conformal  mapping  of  Fig[l]  Fg-IU 


Fig.  3.  The  harmonic  skeleton  for  Fig [T] 


Fig.  8.  Area-preserving  flattening  of  Fig0 
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6  Conclusions 

In  this  paper,  we  presented  a  new  method  for  flattening  of  branched  medical 
surfaces  based  on  the  theory  of  optimal  mass  transportation.  Starting  from  a 
conformal  equivalence,  we  constructed  a  density  map  and  then  generated  an  area¬ 
preserving  representation  of  the  original  3D  surface.  Potential  applications  in¬ 
clude  the  detection  and  visualization  of  pathologies  such  as  stenoses  and  polyps. 
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